COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	\magnification\magstephalf
C00038 ENDMK
C⊗;
\magnification\magstephalf
\input macro.tex
\def\today{\ifcase\month\or
  January\or February\or March\or April\or May\or June\or
  July\or August\or September\or October\or November\or December\fi
  \space\number\day, \number\year}
\baselineskip 14pt
\def\fnc#1{\mathop{\rm #1}\nolimits}
\rm
\line{\sevenrm a01.tex[106,phy] \today\hfill}

\bigskip
\line{\bf Case Study: A Historical Computation of Pi\hfil}
\line{\it ``Let's have another cup of coffee, and let's have another piece
of $\pi$''\hfil}

\smallskip
In 1706 the English mathematician  Machin (1680-1752) published the  value
of $π$ computed  to 100  decimal places.  Previous  calculations had  mainly
been based  on approximating  the circle  by polygons  with a  very  large
number of  sides;  the  best of  these  was  that of  Ludolph  van  Ceulen
(1540--1610), correct to 35 decimal places.

Machin may have hoped to discover a repeating pattern in the digits of $π$.
We now know that $\pi$ is irrational, so that  the digits don't repeat.   For
most practical applications, ten digits or so of $\pi$ exceed the accuracy  of
any  current  methods  of  measurement.   If  we  calculate  the   earth's
circumference using a value of $\pi$ correctly rounded to ten decimal  places,
the maximum error  in the result  is $8000 \hbox{ miles} \times 0.5\times 10↑{-10}$, 
which  is less than  a millimeter  (0.026~in.). 
According to  Pipkin and  Ritter, {\it Science\/} Vol.~219, 
No.~4587, 25~Feb.~1983, 
the most accurate  measurement involving length is probably of the speed of light, 
to about 8~digit accuracy.

\vfill\eject

Machin's method was based on a division of the circle.  First cut the unit
circle into eight pie-shaped slices (sectors).  Each has a portion of  the
circle's circumference  of  length $\pi/4$,  and  has a  central  angle  whose
tangent is~1.   If  you can  calculate  the  inverse tangent  of~1  very
accurately, you can multiply the result by~4 and get~$\pi$.

Not long before  (1670),  James Gregory  had found a power  series for the 
inverse tangent function.  It was 
$$\arctan (x) = x-x↑3/3 + x↑5/5 - x↑7/7 + x↑9/9\ldots\;,\hbox{ etc.}$$
Now that calculus has been invented, we can get Gregory's  formula easily 
by integrating both sides of the equation
$$1/(1 + x↑2) = 1 - x↑2 + x↑4 - x↑6 + x↑8\ldots\;.$$
Gregory's series, however, is divergent if  $x > 1$, and for  $x = 1$ is very 
slowly  convergent.  If we  calculated a billion  terms of the  series per 
second, in a hundred years we would have $\arctan (1)$ correct to 15~decimal 
places.  It would take $10↑{85}$ times as long to get it to 100~decimal places 
that way.

Machin instead subdivided the sector, as shown below.  He found
five angles which added up  to $\pi/4$.  Four of them  were equal, and all  of
them had simple numbers as tangents: $1/5$, $1/5$, $1/5$, $1/5$, 
and $-1/239$.   You
can confirm by standard trigonometric  identities for sums of angles  that
the tangent of the sum is~1, or evaluate the complex number $(5+i)↑4(239-i)$.
To get $\pi/4$, then, we need only  evaluate $4\arctan (1/5) - \arctan (1/239)$,
using Gregory's series.

To get $\pi$ to 100 places, the error in $\pi$ must be less than  
$0.5\times 10↑{-100}$; the error in  $\pi/4$ must be less than  
$0.125\times 10↑{-100}$, and  the error  in $\arctan (1/5)$ must be less than  
$0.03125\times 10↑{-100}$.  To do this, we must keep the  effects of  about 
100~rounding errors  from adding  up to  
$3\times 10↑{-102}$, so working with numbers calculated to  105~decimal places  
should suffice.  Also, we must use enough terms in the sequence
$${1\over 5}\;,\;{1\over 3\times 5↑3}\;,\;
{1\over 5\times 5↑5}\;,\;{1\over 7\times 5↑7}\;,\;
{1\over 9\times 5↑9}\;\cdots$$
that the next one is less than   $3\times 10↑{-102}$;
it is sufficient to stop at $1/(141\times 5↑{141})$, the 71-st term.  
Similarly,  the calculation of 
$\fnc{arctan}(1/239)$ may stop with $1/(41\times 239↑{41})$, the 21-st term.

To show how Machin must have done  the calculation, here we work out $\pi$ to
10~decimal places, using intermediate numbers calculated to twelve places.

\vfill\eject

\noindent
$$\halign{$\hfil #$&${}#\hfil $&$\hfil #$&${}#$&$\quad #$\qquad%
&$\displaystyle{#}$\hfil\cr
\noalign{\hbox{(1) Calculate $\arctan\, (1/5)$}}
\noalign{\medskip}
1 & \div 5 = 0.20000\ 00000\ 00 & \div 1 & = 0.20000\ 00000\ 00 & (+)
&\left({1\over 5}\right)\cr
&&\hbox{current sum}& = 0.20000\ 00000\ 00&\cr
  & \div 5 = 0.04000\ 00000\ 00&\cr
  & \div 5 = 0.00800\ 00000\ 00	& \div 3 & = 0.00266\ 66666\ 66	& (-)
&\left({1\over 3\cdot 5↑3}\right)\cr
&&\hbox{current sum}& = 0.19733\ 33333\ 34\cr
  & \div 5 = 0.00160\ 00000\ 00\cr
  & \div 5 = 0.00032\ 00000\ 00 & \div 5 & = 0.00006\ 40000\ 00	& (+)
&\left({1\over 5\cdot 5↑5}\right)\cr
&&\hbox{current sum}& = 0.19739\ 73333\ 34\cr
  & \div 5 = 0.00006\ 40000\ 00\cr
  & \div 5 = 0.00001\ 28000\ 00 & \div 7 & = 0.00000\ 18285\ 71 & (-)
&\left({1\over 7\cdot 5↑7}\right)\cr
&&\hbox{current sum}& = 0.19739\ 55047\ 63\cr
  & \div 5 = 0.00000\ 25600\ 00\cr
  & \div 5 = 0.00000\ 05120\ 00 & \div 9 & = 0.00000\ 00568\ 88 & (+)\cr
&&\hbox{current sum}& = 0.19739\ 55616\ 51\cr
  & \div 5 = 0.00000\ 01024\ 00\cr
  & \div 5 = 0.00000\ 00204\ 80 & \div 11 & = 0.00000\ 00018\ 61 & (-)\cr
&&\hbox{current sum}& = 0.19739\ 55597\ 90\cr
  & \div 5 = 0.00000\ 00040\ 96\cr
  & \div 5 = 0.00000\ 00008\ 19 & \div 13 & = 0.00000\ 00000\ 63 & (+)\cr
&&\hbox{current sum}& = 0.19739\ 55598\ 53\cr
  & \div 5 = 0.00000\ 00001\ 63\cr
  & \div 5 = 0.00000\ 00000\ 32 & \div 15 & = 0.00000\ 00000\ 02 & (-)\cr
&&\hbox{current sum}& = 0.19739\ 55598\ 51\cr
\noalign{\bigskip}
\noalign{\hbox{(2) Calculate $\arctan\, (1/239)$}}
\noalign{\medskip}
1 & \div 239 = 0.00418\ 41004\ 18 & \div 1 & = 0.00418\ 41004\ 18 & (+)\cr
&&\hbox{current sum}& = 0.00418\ 41004\ 18\cr
  & \div 239 = 0.00001\ 75066\ 96\cr
  & \div 239 = 0.00000\ 00732\ 49 & \div 3 & = 0.00000\ 00244\ 16 & (-)\cr
&&\hbox{current sum}& = 0.00418\ 40760\ 02\cr
  & \div 239 = 0.00000\	00003\ 06\cr
  & \div 239 = 0.00000\ 00000\ 01 & \div 5 & = 0.00000\ 00000\ 00 & (+)\cr
&&\hbox{current sum}& = 0.00418\ 40760\ 02\cr}$$

Machin's algorithm, even though it  works with 100-digit numbers, is  kept
reasonable because  it avoids  fullscale multiplications and divisions  of
long numbers  by each  other.  The  operations used  on long  numbers  are
addition, subtraction, multiplication  by~4, and  division by integers  no
greater than~239.  Even so, Machin, working by hand, must have  calculated
about 27000 digits  of quotients.   Notice that, if  he had  tried to  get
$(1/239)↑2$ by multiplying  $1/239$ by itself, that alone  would have involved
10000 digits, and to  do the whole computation  that way would have  taken
years.
	
Let's look in more detail at Machin's  calculation of the powers of $1/239$.
To do it  by hand,  the methods  are obvious.  What  if we  have a  pocket
calculator capable of  working with  ten digit  numbers only;  how can  we
organize the calculation  to use it?   Let's jump into  the middle of  the
calculation, to 30  decimal places;  we already  have $1/239$,  and we  want
$(1/239)↑2$.
$$1/239=\,\hbox{0.00418 41004 18410 04184 10041 84100}.$$
Let's  scale  the  numbers  up by a  factor of  $10↑{30}$,  so we're working 
with integers.  We want  to divide  00418  41004  18410  04184  10041  84100 
by~239.   We'll  develop  the quotient  in  5~digit  groups
(bigits using base 100000),  starting at  
the left. Dividing 418 by 239, we get a quotient of~1, and a remainder of
$418 - 1\times 239 = 179$. The remainder, 
because of its position in the number,  represents the value 
$179\times 10↑{25} = 179 00000\times 10↑{20}$.  The next group  of digits 
in the dividend, 41004, represents  $41004\times 10↑{20}$; putting the 
two together by the formula
$$179\times 10↑5  + 41004 = 17941004\,,$$
we divide this number by  239  to get as quotient the next five  digits of
$1/239$: 75066, with remainder~230.  Continuing this way,  we develop  the
quotient as shown below, where no calculation step involves more than  ten
digit numbers,  and  each  step  is  done  exactly,  without  rounding  or
discarding remainders.

$$\vbox{\tabskip=0pt\offinterlineskip
\halign{\strut$\lft{#}$&#&#&\hfill#\quad
&\hfill#\quad&\hfill#\quad&\hfill#\quad&\hfill#\quad
&\hfill#\quad&\hfill#\quad&\hfill#\quad&\hfill#\cr
&&&1&75066&96311&33908&72008&54326\cr
\omit&&&\multispan6\hrulefill\cr
&\quad 239\quad&\vrule height12pt depth4pt&\quad 418&41004&18410&04184&10041&84100\cr
\omit&\hrulefill\cr
239\times 1=&&&239\cr
\omit&&\multispan2\quad\hrulefill\quad\cr
\noalign{\smallskip}
&&&179&41004\cr
\noalign{\smallskip}
239\times 75066=&&&179&40774\cr
\omit&&&\multispan2\quad\hrulefill\quad\cr
\noalign{\smallskip}
&&&&230&18410\cr
\noalign{\smallskip}
239\times 96311=&&&&230&18329\cr
\omit&&&&\multispan2\quad\hrulefill\quad\cr
\noalign{\smallskip}
&&&&&81&04184\cr
\noalign{\smallskip}
239\times 33908=&&&&&81&04012\cr
\omit&&&&&\multispan2\quad\hphantom{0}\hrulefill\quad\cr
\noalign{\smallskip}
&&&&&&172&10040\cr
\noalign{\smallskip}
&&&&&&172&09912\cr
\omit&&&&&&\multispan2\quad\hrulefill\quad\cr
\noalign{\smallskip}
&&&&&&&129&84100\cr
\noalign{\smallskip}
&&&&&&&129&83914\cr
\omit&&&&&&&\multispan2\quad\hrulefill\quad\cr
\noalign{\smallskip}
&&&&&&&&186\cr}}$$

\noindent
So the  quotient is  $(1/239)↑2=$ 0.00001 75066  96311, etc.  The  algorithm
embodied in the above calculation can be written down explicitly.
The 5-digit groups of the dividend are in variables 
{\tt DVND}$[1\ldt 6]$.
We will put  the quotient digit  groups into the  result variables  
{\tt RES}$[1\ldt 6]$.
We  will use other  variables:  {\tt PAR} to  hold the  partial
dividend, {\tt QUO} for  the quotient, and  {\tt REM} 
the remainder, of the  calculation steps.

\vfill\eject

The algorithm goes like this:

\smallskip
\disleft 40pt:: 
Set $\hbox{\tt REM}= 0$ ($\ast$ Initializing $\ast$).

\smallskip\noindent
Do the following four steps for each value of ${\rm I} = 1,2,\ldots,6$.

\smallskip
\disleft 40pt::
Set $\hbox{\tt PAR}=\hbox{\tt REM}\times 10↑5 + \hbox{\tt DVND}↓{\rm I}$.

\smallskip
\disleft 40pt:: 
Set $\hbox{\tt QUO}=\hbox{the integer part of {\tt PAR}}/239$.

\smallskip
\disleft 40pt:: 
Set $\hbox{\tt REM}=\hbox{\tt PAR}- 239 \times\hbox{\tt QUO}$.

\smallskip
\disleft 40pt:: 
Set $\hbox{\tt RES}↓{\rm I} =\hbox{\tt QUO}$.

In the above algorithm, we use at each step a remainder from the  previous
step.  Instead of having  separate rules for the  first step (${\rm I}=1$), it  is
simpler to say that there is a remainder, but it is zero.

If Machin did his calculation on  paper, it could have appeared much  like
the one above.  If  he used a blackboard  with limited space, however,  he
could have calculated the arctans without keeping all the numbers for  the
whole time.   At the stage of the calculation that adds $1/(N\times 239↑N)$ 
to the  current  sum,  all  the (human) calculator  needs is the  values of~$N$,
$1/239↑N$, $1/(N\times 239↑N)$, and the current sum; 
all other  numbers can be  erased.  To go 
on to the next  value of $N$, the calculator does this:

\smallskip
\disleft 40pt:: 
Increase $N$ by 2.

\smallskip
\disleft 40pt:: 
Divide the old value of $1/239↑N$  by 239, twice; 
	the result is the new value.

\smallskip
\disleft 40pt:: 
Divide $1/239↑N$  by $N$ to get  $1/(N\times 239↑N)$.

\smallskip
\disleft 40pt:: 
Increase or decrease the current sum by  $1/(N\times 239↑N)$.

\smallskip
\noindent
Using symbolic execution:
$$\vcenter{\halign{#\hfill\quad&\hfil$\displaystyle{#}$\hfil\quad%
&\hfil$\displaystyle{#}$\hfil\quad&#\hfill\cr
\hfill COUNT&\hbox{POWER}&\hbox{TERM}&\hfil SUM\cr
\noalign{\medskip}
\hfill $N$&{1\over 239↑N}&{1\over N\times 239↑N}&Sum of first $N$ terms\cr
\noalign{\smallskip}
Increase $N$ by 2\cr
\noalign{\medskip}
\hfill $N+2$\cr
\noalign{\smallskip}
Divide POWER by 239\cr
&{1\over 239↑{N+1}}\cr
\noalign{\smallskip}
Divide POWER by 239 again\cr
&{1\over 239↑{N+2}}\cr
\noalign{\smallskip}
Set TERM $=$ POWER/COUNT\cr
&&{1\over (N+2)\times 239↑{N+2}}\cr
\noalign{\smallskip}
Set SUM $=$ SUM $≠$ TERM&&&Sum of first $N+2$ terms\cr}}$$

\vfill\eject

In recent years, calculations of $\pi$ by computer have reached astronomical
scales.  A report by D.~Shanks and J.~W. Wrench, in {\sl The Mathematics of
Computation} vol.~16 (1962), 76--99, contains the first 100,000 digits.  
In~1967, Jean Guillond computed $\pi$ to 500,000 places. [RWF: reference]

In 1984, mathematicians using a very large and fast computer at the
University of Tokyo computed pi to 16~million decimal places, a~number
so large it would fill two Sunday editions of a major newpaper.
[Reference: {\sl The Mathematics of Computation\/}, 1984?]

\bigskip

\disleft 38pt::
{\bf Exercises:}

\smallskip
\disleft 38pt:(1):
Compute $e$ ($=2.71828\ldots$) to 100 decimal places.

\smallskip
\disleft 38pt:(2):
Compute $1/\pi$ to $100{\rm D}$.

\smallskip
\disleft 38pt:(3):
Compute $\sqrt{\pi}$ to $100{\rm D}$.  (Hint: divide 
$\pi$ by numbers like $1.1↑2$, $1.01↑2$, $1.001↑2$,
until it reaches $1.000\ldots$; keep the product of the square roots 
of the divisors.)

\smallskip
\disleft 38pt:(4):
Compute $\pi$ by the formula $6\arccos(1/2)$.

\smallskip
\disleft 38pt:(5):
Compute $e↑{\pi}$ to $100{\rm D}$. (Quite hard.)

\smallskip
\disleft 38pt:(6):
Generalize the procedures used in the $\pi$ computations, for inputs from arbitrary
locations.

\smallskip
\disleft 38pt:(7):
Compute the natural logarithms of the first five prime numbers, 
2, 3, 5, 7 and~11, to
100 decimal places.  Natural logarithms of numbers between 0 and~2 can be
calculated by the expression
$$\ln(1+x) = x-x↑2/2+x↑3/3-x↑4/4+x↑5/5\cdots\,.$$
If we can express each of the five primes as a product of several numbers close
to~1, we can add up the logarithms of these numbers to get the desired results.
We could use, say 
$$2 = 3/2 \times 4/3\,,\hbox{ so }\ln 2 =\ln(3/2)+\ln(4/3)\,.$$

\smallskip
\disleft 38pt::
A better choice of numbers close to 1 is:
$$\eqalign{A&= 1+1/32  = 33/32\cr
B&= 1+1/63  = 64/63\cr
C&= 1+1/27  = 28/27\cr
D&= 1+1/80  = 81/80\cr
E&= 1+1/242 = 243/242\cr
\noalign{\vskip 4pt}
2&= A↑{10}  B↑7  C↑7  E↑5\cr
3&= A↑{16}  B↑{11}  C↑{11}   E↑8\cr
5&= A↑{24}  B↑{16}  C↑{16}  D↑{-1}  E↑{12}\cr
7&= A↑{28}  B↑{19}  C↑{20}  E↑{14}\cr
11&= A↑{35}  B↑{24}  C↑{24}  E↑{17}\cr}$$
Then $\ln 2 = 10\ln A + 7\ln B + 7\ln C + 5\ln E$, etc.

\smallskip
\disleft 38pt::
An alternate method calculates the natural logarithms of $1/A$, $1/B$, etc;
$(1/A = 1-1/33)$ by the series expansion
$$\ln(1-x) = -(x+x↑2/2+x↑3/3+x↑4/4\cdots)\,,$$
avoiding the alternating signs.

\smallskip
\disleft 38pt::
To check your results, the 91st to 100th decimal places of $\ln 2$ are
[RWF: fill in]


\bigskip
\line{\copyright 1984 Robert W. Floyd; 
First draft (not published) February 23, 1984\hfil}
%revised: Date; subsequently revised.\hfill}

\bye